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In this paper we make an effort to understand the interaction of turbulence generated by the magnetorotational instability 
(MRI) with turbulence from other sources, such as supernova explosions (SNe) in galactic disks. First we perform a 
linear stability analysis (LSA) of non-ideal MRI to derive the limiting value of Ohmic diffusion that is needed to inhibit 
the growth of the instability for different types of rotation laws. With the help of a simple analytical expression derived 
under first-order smoothing approximation (FOSA), an estimate of the limiting turbulence level and hence the turbulent 
diffusion needed to damp the MRI is derived. Secondly, we perform numerical simulations in local cubes of isothermal 
nonstratified gas with external forcing of varying strength to see whether the linear result holds for more complex systems. 
Purely hydrodynamic calculations with forcing, rotation and shear are made for reference purposes, and as expected, 
non-zero Reynolds stresses are found. In the magnetohydrodynamic calculations, therefore, the total stresses generated 
are a sum of the forcing and MRI contributions. To separate these contributions, we perform reference runs with MRI- 
stable shear profiles (angular velocity increasing outwards), which suggest that the MRI-generated stresses indeed become 
strongly suppressed as function of the forcing. The Maxwell to Reynolds stress ratio is observed to decrease by an order 
of magnitude as the turbulence level due to external forcing exceeds the predicted limiting value, which we interpret as a 
sign of MRI suppression. Finally, we apply these results to estimate the limiting radius inside of which the SN activity can 
suppress the MRI, arriving at a value of 14 kpc. 
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1 Introduction 

Most differentially rotating astrophysical disks are suscepti- 
ble to the magnetorotational instability (MRI). In accretion 
disks this instability is likely to be the most important source 
of turbulence and angular momentum transport leading to 
mass accretion (originally proposed by Balbus & Hawley 
1991). It has also been proposed that the MRI would play a 
role in driving turbulence in the outer parts of galactic disks 
where the stellar forcing is of minor importance (Sellwood 
& Balbus 1999; see also Piontek & Ostriker 2007 and refer- 
ences therein). 

Whether the powerful stellar forcing suppresses or in- 
teracts with the galactic MRI in the inner parts of galac- 
tic disks is still an open question. The study by Workman 
& Armitage (2008), who investigated the interaction of the 
MRI with nonhelically forced turbulence in isothermal non- 
stratified local cubes of gas with Keplerian differential rota- 
tion, can be regarded as a first step to understand this issue. 
They found that external turbulence can indeed modify the 
turbulent transport when the level of turbulence becomes 
greater than the strength of the MRI-generated turbulence, 
and that the degree of the influence is dependent on the 
forcing scale. Only turbulence forced at large wavenumbers 
(smallest scales) was interpreted to suppress MRI, while 
large-scale forcing seemed to enhance it. In this study we 
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show that the forcing wavenumber dependence can be un- 
derstood with the help of the linear stability analysis. We 
also show that reference calculations with MRJ-stable rota- 
tion profiles are needed to fully determine the level of tur- 
bulence and angular momentum transport that results from 
the MRI in comparison to other sources of turbulence. 

The paper is organised as follows. In Sect. |2] we present 
a linear stability analysis of non-ideal MRI to derive the 
limiting value of Ohmic diffusion that is needed to inhibit 
the growth of the instability for different types of rotation 
laws. We also present an estimate of the level of turbulence 
needed to produce turbulent Ohmic diffusion strong enough 
to suppress MRI in the isotropic homogeneous case using 
the first-order smoothing approximation (hereafter FOSA). 
In Sect. [3] we introduce the model that we use to study the 
dynamics of the MRI under the influence of external non- 
helical forcing. First we apply this model to carry out sim- 
ple one-dimensional calculations to reproduce the linear re- 
sults of Sect. |2] and then to three-dimensional hydro- and 
magnetohydrodynamic models of isothermal, nonstratified, 
rotating and shearing cubic domains subject to rotation laws 
both unstable and stable with respect to the MRI. The latter 
have been used to calibrate the former to get a better under- 
standing of turbulence resulting from the MRI with respect 
to the turbulence resulting from the interaction of external 
forcing, rotation and shear. The results of these calculations 
are presented in Sect.|4] In Sect.|5]we present a simple appli- 
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cation of our results to galactic disks. In Sect. |6]we present 
our conclusions. 



From these equations the following dispersion relation is 
obtained: 



2 Linear stability analysis of non-ideal MRI 

Jin (1996) presented linear stability analysis of non-ideal 
MRI for disks with Keplerian rotation profiles. Here we 
generalise the theory to include any rotation profile of the 
form il oc r^"^, further assuming that the system is incom- 
pressible and pressureless. 

We choose a local Cartesian frame of reference, in 
which X, y and z denote the radial, toroidal and vertical di- 
rections, respectively. We start by linearizing the pressure- 
less ideal equation of motion 

^ + (U- V)U + 2nxU --J X B = , 

ot p 

and the non-ideal induction equation 



--VxiUxB) 



(1) 



(2) 



around the location of the centre of the box, ro, in a differen- 
tially rotating disk that is subject to a weak uniform vertical 
magnetic field BqEz. Here U is the velocity field, B is the 
magnetic field, J = /Iq ^ V x B is the current density, /io is 
the vacuum permeability, and rj is the Ohmic diffusion co- 
efficient. The angular velocity at ro is flo, and it varies with 
the distance from the centre as cx r^''. Such differen- 
tial rotation translates into a radially linear shearing veloc- 
ity Uq — —qfloxey after the linearization of Eq. ([T]i- Small 
fluctuations u = [ua:,Uy,Uz] , b = [bx,by,bz] are now 
added to the equilibrium state Uq = [0, —qilox, 0] , Bq = 
[0, 0, Bq]. We consider only fluctuations in the vertical (z) 
direction, as fluctuations of this type have been shown to 
have the highest growth rates (Balbus & Hawley 1992); in 
this case the equations for the z-components of velocity and 
magnetic field decouple from the horizontal ones and there- 
fore they can be neglected. The set of linearised equations 
for the horizontal fluctuating fields in component form read 

dux „ „ Bq dbx 

— — = 2ilQUy H — , 

ot fioP oz 

duy Bq dby 

ot hqp az 



(3) 



(4) 



dt 
db 



B, 



dux 
oz 



d^bx 



dz'^ 
du-u 



d'b. 



(5) 



(6) 



ot Oz Oz'^ 

Let us now assume that the solutions to these equations are 
of the form /, = /j-e*'-'^^""*''. Substitution yields a set of 
algebraic equations 



lUJUx 



-iubx 
-iub,, 



2 floiiy 



■1 ^0 u 

ik Ox , 

fJ-oP 



Bn 



= — (2 — g) floiix + ik by , 



ikBoUx - k rjbx , 
—qVlobx + ikB^Uy — k^r/by 



(7) 

(8) 

(9) 
(10) 



- (2 (2 - q) nl + + 2k^vl) uj^ 

~ in^ (4 (2 -q)nl + 2k'^vD lu 

+ k'^vl{k'^vl~2qnl) + 2{2-q) ^1^1^ = 0, 



(11) 



where vj. — and i}„ 

A Mop " 



r/k^. We have numerically 



solved the dispersion relation searching especially for the 
limiting Ohmic diffusion rate at which all the MRI modes 
become damped. In practise we solve for all the four roots 
of the dispersion relation and require that all the imaginary 
parts are negative. The results of this analysis are plotted in 
Figs.[I]and|2l 

In Fig. [T] we show all imaginary parts of the four roots 
of the dispersion relation for q = 1.5, kj^ax/k = 1, where 

fcmax = 2 '\fq~~^' functions of the Ohmic diffusion 
rate ilm. We have scaled ilm with the maximum growth 
rate, Fmax, of the ideal MRI. From these results it can be 
seen that in the ideal limit = the growth rate of the 
MRI mode is positive and equal to Fmax = |fio- As the dif- 
fusion is increased, the growth rate decreases and reaches 
zero at Jim w l.OFmax- As the growth rates of the other 
modes are negative for all values of ft^, then this is the 
limiting diffusion rate above which the instability will be 
damped for k = fcmax- Our results are in agreement with 
those of Jin (1996); due to the incompressibility assump- 
tion the inertial wave mode is absent in our analysis, but the 
limiting Ohmic diffusion rate for MRI damping is identical 
for the Keplerian case in both studies. 

In Fig. 12] we have determined the limiting Ohmic diffu- 
sion rate for each value of q, < q < 2, where the system 
is ideally unstable, for all values of A: = [0, fccrit], where 
fccrit = ^ \/2q is the largest MRI-unstable ideal wavenum- 
ber. Analytical 'stability criterion' can also be derived from 
the analysis: MRI will be damped if 



kVA 



2{2-q) 



(12) 



At the limit q ^ 2, the Ohmic diffusion rate becomes in- 
finite, meaning that in this regime the flow cannot be sta- 
bilized by the Ohmic diffusion. We find that for small q the 
results for the maximally excited and all unstable wavenum- 
bers are identical, whereas for larger q they start to devi- 
ate more and more indicating that some wavenumbers are 
harder to damp than the one with the largest growth rate. 
From this plot we can also see that for very small q the lim- 
iting value approaches the maximum growth rate, Fmax, of 
the ideal MRI. The limiting Ohmic diffusion rate increases 
with q above Fmax so that for rotation profiles near the 
Rayleigh unstable regime (q > 2) the Ohmic diffusion rate 
needed to damp the MRI is roughly three times Fmax- 
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Fig. 1 Imaginary parts of the roots of the non-ideal disper- 
sion relation for the Keplerian case q — 1.5 as function of 
the Ohmic diffusion rate scaled with the maximum growth 
rate Fmax of the MRI. Solid thick line represents the (fast) 
MRI mode, dashed line the decay mode, and the dotted line 
the slow mode (given by two roots with equal imaginary 
parts). The dashed-dotted line indicates ijjj/i7o=0. 



Fig. 2 Limiting Ohmic diffusion rate fim at which 
MRI modes corresponding to the range of wavenumbers 
k = [0, /cciit] (solid line with diamonds) and for a single 
wavenumber k = fc,„ax (dashed line with triangles) are 
damped as functions of q. The results are scaled with the 
maximum growth rate of the ideal MRI, Fuiax- The dotted 
line shows the limit at which the Ohmic diffusion rate and 
maximum growth rate of MRI are equal for reference. 



2.1 FOSA-prediction for isotropic homogeneous 

turbulence 3 Model 



Let us now use the information obtained from the linear sta- 
bility analysis to make a prediction for the nonlinear setup 
we use in the rest of the paper. Firstly we need an estimate 
of the turbulent Ohmic diffusion resulting from the isotropic 
external forcing at a wavenumber kf. In order to obtain this 
we can use the first-order smoothing approximation (FOS A) 
result 



Vt 



3 ' c "Tins 



(13) 



where we have assumed that the Strouhal number, St = 
TcUimsfcf , is equal to unity. This value has been confirmed in 
the kinematic regime numerically by Sur et al. (2008) from 
non-shearing setups using the test field method (Schrinner et 
al. 2005, 2007). In the presence of shear the value of rjt can 
increase slightly but is still of the same order of magnitude 
(Brandenburg et al. 2008; Mitra et al. 2009). Therefore we 
obtain 



(14) 



for the turbulent Ohmic diffusion rate. From the linear sta 
bility analysis we know that MRI will be damped if 

qAiq) 



max " 



(15) 



where A{q) varies from 1 to roughly 3 as q is increased from 
to 2 (see Fig. |2]l. By substituting the FOSA-expression, 
Eq. (fTSl l. we can simplify this further, yielding the condition 



> 



3qA{q)n 
2kf ' 



(16) 



for the suppression of MRI. Note that for a given Umis this 
condition is dependent on the spatial scale of the turbulence 
so that the larger the scale, the harder it is to damp MRI. 



We use a model setup similar to that used in Liljestrom et al. 
(2009), except that here we add an external forcing function 
in the Navier-Stokes equation. The computational domain 
is a cube of size (27r)'^, filled with isothermal gas; verti- 
cal gravity and stratification are neglected. The calculations 
are local and utilise the shearing-sheet approximation. The 
equations to be solved read 



Vp 

vt 

Vu 



— = -(it • V)m - qfloUxey - 2 fioe^ x u 



--Vp- 
P 



-J xB 

P 



J force 1 



■DA 

'vt 



= u X B + qrioAyBx + 77hV^A 



(17) 



(18) 



(19) 



where V/Vt = d/dt + Uod/dy includes advection by the 
shear flow, u is the departure from the mean flow Uq, p 
is the density, p is the pressure, A is the magnetic vec- 
tor potential, B = V x A is the magnetic field, and 
J = V X B//io is the current density. In order to get as 
close as possible to the ideal limit, we use hyperviscous op- 
erators of the form to replace the ordinary viscosity V'^- 
operators; i/h and stand for the hyperviscous kinematic 
viscosity and magnetic diffusivity, respectively. With this 
procedure we aim at maximizing the Reynolds number in 
the quiescent regions of the flow while diffusing and damp- 
ing fluctuations near the grid scale. Compared to direct sim- 
ulations (e.g. Haugen & Brandenburg 2006) with uniform 
viscosities, smaller grid size can be used to resolve the flow. 
The forcing function /force is given by 



fix, t) = Re{iV/fc(t) exp[ifc(i) • x - 



(20) 
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where x is the position vector, N — foCs{kcs/St)^^^ is a 
normalization factor, /o is the forcing ampHtude, k ~ \k\, 
St is the length of the time step, and — tt < <j){t) < tt a 
random delta-correlated phase. The vector fk is given by 

fc X e 

= 7P=triF ' ^^^^ 

where e is an arbitrary unit vector. Thus, fk describes non- 
helical transversal waves with = 1, where k is chosen 
randomly from a predefined range in the vicinity of the av- 
erage nondimensional wavenumber k[ /ki at each time step. 
Here ki is the wavenumber corresponding to the domain 
size, and kf is the wavenumber of the energy-carrying scale. 

Periodic boundary conditions are applied in all three di- 
rections; in the radial direction we account for the shear flow 
Uq by making use of the shearing-sheet approximation (e.g. 
Wisdom & Tremaine 1988): 

/(iL,, y, z) = /(-iL,, y + qn^L^t, z), (22) 

where / stands for any of the seven independent variables, 
Lx for the radial extent of the computational domain, and t 
is the time. 

The domain is initially threaded by a weak magnetic 
field. 



A ~ AoBy COS kAX cos kAy cos kAZ 



(23) 



Thus, the magnetic field contains periodic x- and z- 
components with amplitude Aq. The values of /ca, and 
Aq are selected so that both the wavenumber with the largest 
growth rate, fcmax, and the largest unstable wavenum- 
ber, feci it, are well resolved by the grid; in practise this 
means that in the three-dimensional calculations we adopt 
kA/ki = 1, flo — 0.2csfci and Aq ~ 0.2y/JIopoCsk^^, re- 
sulting in kniax = 0{ki). For the initial setup, the other 
condition for the onset of MRI, namely /? ^ 1, where 
(3 — 2fj,op/BQ is the ratio of the thermal to magnetic pres- 
sure, is also satisfied as /? is at least 50 at the maximum 
values of the magnetic field. 

For all the simulations we use the PENCIL COD^H 
which is a high-order (sixth order in space, third order in 
time), finite-difference code for solving the MHD equations 
(e.g. Brandenburg & Dobler 2002). 

The local calculations have been carried out at two dif- 
ferent resolutions, namely 128 (in ID) and 64'^ (in 3D). In 
the three-dimensional calculations the hyperdiffusion coef- 
ficients used are i/h = % = 2 x 10~'^ Cgk^^ . 

4 Results 

4.1 Linear regime 

To validate the numerical method used in the three- 
dimensional nonlinear calculations, we first try to reproduce 
the linear results of Sect. |2] with a one-dimensional setup 
using ordinary magnetic diffusion with -operator and a 
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Fig. 3 The growth rate uj of the MRI-mode with fcmax = 
fci, rio/icski) = 1, and va/cs = 1 as a function of il^n- 
Solid line: linear stability result for the Keplerian rotation 
law with q = 3/2. Dashed line: the same for flat galactic ro- 
tation curves with q — 1. Triangles and diamonds show the 
growth rates measured from the one-dimensional numerical 
model. 

constant coefficient ry instead of the hyperdiffusion scheme. 
For these runs we choose kmax = ki, ilo/{cski) — 1, 
va/cs — 1, vary q and 77, and impose a uniform magnetic 

field of strength Bo/{,/p^Cs 



_ 2 
~ 2 



1 in the verti- 



cal (z) direction at each time step. This procedure forces 
the system to develop MRI at the wavenumber k = fcmax 
with the growth rate Fmax at the ideal limit (fim = 0). In 
Fig. [3] we plot the growth rate of the MRI mode as function 
of rim/rmax for two valucs of q. As can be seen from this 
figure, the results of both methods are in perfect agreement. 

We also compute the time averages of the Maxwell and 
Reynolds stresses, A/^jy — (J-q^ i^xby) and R^y — {u^Uy), 
from these calculations. In the limit il„i = the linear stress 
ratio reads 



M. 



xy 



A-q 



Rx 



(24) 



http://www.nordita.org/software/pencil-code/ 



1^2; y q 

as derived by Pessah et al. (2006). As can be seen from 
Fig. |4] where we plot the stress ratio as a function of Vl^ 
for the Keplerian {q = 3/2) and flat (q = 1) rotation laws 
from the one-dimensional calculations, our results agree 
with Eq. (|24] | for = 0. As we increase rj and thereby 
r^m, the stress ratios decrease, reaching minima somewhat 
before the limiting Ohmic diffusion rates above which the 
MRI is damped and the growth rates are negative. Although 
the stress ratios in the fully nonlinear state of the MRI differ 
from the linear ones (see e.g. Liljestrom et al. 2009), these 
results suggest that the sress ratio should show a decreasing 
trend as a function of the turbulent Ohmic diffusion rate also 
in the nonhnear regime. 

4.2 Three-dimensional calculations with external 
forcing 

Let us now consider the implications of the above analy- 
sis to the earlier numerical work done by Workman & Ar- 
mitage (2008). They used a setup very similar to the one 
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Fig. 4 The Maxwell to Reynolds stress ratio —Mxy/Rxy 
as a function of fim obtained from the one-dimensional cal- 
culations with Keplerian rotation law (q = 3/2, solid line) 
and with flat rotation curve (q = 1, dashed line). 

presented in Sect. [3] varying the amplitude /o and the en- 
ergy carrying wavenumber k{ of the forcing function given 
by Eq. (|20] i. Two different forcing scales were investigated, 
namely a large-scale forcing with fcf/fci = 1.5 and a small- 
scale forcing with fcf/fci = 40, the other relevant param- 
eters chosen to be rio/(csA;i) = 0.2 and q = 1.5. For 
these parameters our analysis predicts roughly a factor of 
30 difference in the limiting rms-values, the exact num- 
bers reading Wrms/cs = 0.6 for the large-scale forcing, and 
"rms/cs = 0.0225 for the small-scale one. Workman & Ar- 
mitage (2008) report a running time average of the kinetic 
energy density E'kin = 0.0012 for the pure MRI case, and 
roughly 15 times more kinetic energy for the largest forcing 
amplitude investigated. Assuming that the average density 
is roughly one, the rms-velocities would read Mrms/cs ~ 
0.049 and 0. 190, for the pure and maximally forced runs, re- 
spectively. Based on this estimate, all of the cases with the 
large-scale forcing investigated by Workman & Armitage 
(2008) produced turbulence whose level was below the sup- 
pression limit; we propose that this is the reason why no 
clear MRI damping was seen with large-scale forcing. On 
the other hand, considering the small-scale forcing, the lim- 
iting rms-velocity is smaller than the one of pure MRI, so 
MRI should be suppressed at any forcing amplitude. 

We perform four sets of three-dimensional nonlinear 
calculations differing by the effective forcing wavenum- 
ber in the range fcf/fci = [1.5,3,5,10]. The parameters 
common to all runs are rio/(csfci) — 0.2, q — 1.5, and 
A{q) « 2. The predicted limiting rms-velocities for the cho- 
sen wavenumbers are uum/cs = [0.6, 0.3, 0.18, 0.09]. In the 
following we calibrate the forcing amplitudes of the calcu- 
lations to produce turbulent rms-velocities below and above 
these limiting values. 

4.2.1 Hydrodynamic calculations 

In a related study (Snellman et al. 2009), the Reynolds 
stresses in forced turbulence with rotation and shear were 
investigated; a non-zero i?^:!, -component naturafly arises in 



Table 1 Hydrodynamic runs with 64"^ mesh points, with 
averages calculated over 150 rotations. Rotation profile with 
(7 = 3/2 i.e. angular velocity decreasing outwards was used 
in these calculations. The number in the brackets in the third 
column indicates the power on ten. 







fo 


{UxUy) /Cg 




fef/fcl 


= 1.5 


0.01 


-1.91 (-6) 


0.07 






0.025 


3.35 (-4) 


0.14 






0.05 


2.25 (-3) 


0.25 






0.10 


7.90 (-3) 


0.40 






0.15 


1.36 (-2) 


0.52 






0.20 


2.07 (-2) 


0.63 






0.25 


2.70 (-2) 


0.72 






0.27 


2.94 (-2) 


0.76 


fcf/fci 


= 3 


0.01 


3.15 (-5) 


0.07 






0.025 


8.62 (-4) 


0.16 






0.05 


3.69 (-3) 


0.27 






0.10 


1.04 (-2) 


0.42 






0.15 


1.70 (-2) 


0.54 






0.20 


2.36 (-2) 


0.64 






0.25 


2.81 (-2) 


0.72 


fcf/fci 


= 5 


0.01 


6.00 (-5) 


0.07 






0.025 


9.98 (-4) 


0.16 






0.05 


3.64 (-3) 


0.26 






0.10 


9.38 (-3) 


0.41 






0.15 


1.41 (-2) 


0.52 






0.20 


1.84 (-2) 


0.61 


fcf/fci 


= 10 


0.01 


6.28 (-5) 


0.06 






0.025 


6.24 (-4) 


0.13 






0.05 


2.22 (-3) 


0.21 






0.10 


5.38 (-3) 


0.34 






0.15 


7.77 (-3) 


0.45 



such a setup. This is in disagreement with the findings of 
Workman & Armitage (2008), who report very small Rxy, 
interpreted to be consistent with zero in their hydrodynamic 
reference runs. This result was taken to show that the hydro- 
dynamic forcing alone cannot produce any stress and that 
all the stresses in the magnetohydrodynamic runs were due 
to the MRI and its interaction with the background turbu- 
lence. This discrepancy between the two very similar stud- 
ies prompted us to perform a series of hydrodynamic runs 
with varying forcing wavenumber and amplitude and mon- 
itor the time evolution of the Reynolds stress component 
Four sets of runs with fcf/fci = [1.5,3,5,10] were per- 
formed. Our results are summarised in Table[T] As is evident 
from this table, we observe a non-zero Rxy in all our cal- 
culations, the magnitude increasing with forcing amplitude 
within each set. No great qualitative differences are found 
between the sets with varying forcing wavenumber These 
results clearly indicate that in the runs with both external 
forcing and MRI, non-zero Reynolds and Maxwell stresses 
can be expected to be generated both from the forcing and 
the MRI, and the two contributions should be carefully sep- 
arated. 

In Snellman et al. (2009) it is shown that for a given 
q, the Reynolds stresses vary as a functions of Q. and S, 
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Fig. 5 Ratio of the Reynolds stresses R^y with positive 
shear parameter to the ones with negative shear parame- 
ters from the hydrodynamic three-dimensional runs as func- 
tion of forcing amplitude /q. Orange color and triangles: 
kf/ki = 1.5, red color and diamonds: kf/ki = 3, green 
color and stars: fcf/fci — 5, and violet color and crosses: 
kf/ki = 10. 

see especially their Fig. 11. In the hydrodynamic run pre- 
sented by Workman & Armitage (2008) the Coriolis num- 
ber, Co — 2Sl/(ui.msfcf), measuring the effect of rotation, 
is roughly 5.3. The results of Snellman et al. (2009) indi- 
cate that in this regime the stress Rxy is small and negative 
although the runs are probably not directly comparable due 
to the different forcing scales. The conclusion of Workman 
& Armitage (2008) that hydrodynamic stresses are small is 
based on this one point, which just by chance lies in an un- 
fortunate part of the parameter space. 

In the study of Snellman et al. (2009) both positive (fl 
decreasing outwards) and negative (fl increasing outwards) 
shear parameters q were investigated systematically. It was 
found that the stresses were not symmetric with respect to 
q in the sense that when changing q —q does not mean 
that the magnitude of the stress remains the same. The same 
behaviour is expected to occur in our setup; to test this, we 
produced sets of calculations identical to the ones presented 
in the previous paragraph but with negative shear parame- 
ter of the same magnitude, namely q ~ —3/2. The results 
are summarised in Table |2l Comparing these numbers to 
the ones listed in Table [T] it can be seen that a significant 
asymmetry indeed exists between the positive and negative 
shear parameter runs, and that this asymmetry is strongly 
dependent on the wavenumber. As can be seen from Fig. |5l 
for the largest wavenumber (kf/ki = 10) investigated, the 
asymmetry is quite weak (the runs with positive shear pa- 
rameter produce Reynolds stresses only slightly larger than 
the negative shear parameter runs) whereas for the smallest 
wavenumber, kf/ki — 1.5, the ratio of positive to negative 
shear parameter-generated Reynols stress is more than two. 

4.2.2 MRI-unstable runs 

Next we turn to the magnetohydrodynamic regime and in- 
vestigate a MRI-unstable setup with k^/ki — \, Vl^ — 



Table 2 The same as Table [T] but rotation profile with 
q = —3/2, i.e. angular velocity increasing outwards, was 
used. 
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Fig. 6 Reynolds stresses R^y generated in the MRI- 
unstable runs with varying forcing wavenumber and forcing 
amplitude, normalised to i?o, the Reynolds stress generated 
in the pure MRI run wihtout forcing. Orange color and tri- 
angles: kf/ki = 1.5, red color and diamonds: fcf/fci = 3, 
green color and stars: fcf/fci = 5, and violet color and 
crosses: fct/^i — 10. The dotted lines plotted over each 
curve indicate the corresponding result from purely hydro- 
dynamic calculations with identical normalisation. 
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Table 3 MRI-unstable runs with nonhelical forcing with resolution 64"^, averages calculated over 250 rotations. The 
mean magnetic field [B) = {{B^y^^ + {By)y^^, where the angular brackets denote time averages and overbar horizontal 
averages. 
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Fig. 7 Maxwell stresses Mxy generated in the MRI- 
unstable runs with varying forcing wavenumber and forcing 
amplitude, normaUsed to Mq, the Maxwell stress generated 
in the pure MRI run wihtout forcing. Colors and symbols as 
in Fig.|6l 



0.2csfci, Aq — 0.2y/fiopo Cgfcf ^, and q = 3/2. For refer- 
ence, we make a pure MRI-run without forcing; the results 
from this run are summarised in the first row of Table |3] 
The typical signatures of MRI (see e.g. Liljestrom et al. 
2009), namely a stress ratio of roughly four and magnetic 
energy exceeding the kinetic energy, are visible. The system 



does not exhibit large-scale dynamo action, so no significant 
mean magnetic fields are present. 

As in the hydrodynamic case, we proceed by mak- 
ing four sets of runs with varying forcing wavenumber 
and amplitude; the parameters for these runs are fcf/fci = 
[1.5, 3, 5, 10] plotted in Figs.lMTTIwith [violet crosses, green 
stars, red diamonds, orange triangles], respectively. Forcing 
amplitudes /o are selected for each set separately so that tur- 
bulence levels below and above the limiting rms-velocities 
for MRI damping uunjcs = [0.6,0.3,0.18,0.09] ai-e ob- 
tained at each fcf ; the values used and the results summarised 
can be found from Table |3] 

In Sect. 14.2.11 it was shown that purely hydrodynamic 
calculations clearly produce non-zero Reynolds stresses 
due to shear, rotation and the adopted forcing function. 
Therefore, in the magnetohydrodynamic regime non-zero 
Maxwell stresses can also be expected to be generated by 
the interaction of a small-scale dynamo and shear without 
the MRI. Therefore, the stresses in a MRI-unstable setup 
can be written as 

Rxy — R^cy ^ R^xy \ (25) 

Mxy = A4f)+M™. (26) 
Here we note that in the study of Workman & Armitage 
(2008) Rxy = Mxy = was apparently assumed in the 
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analysis of the results. To find a definite answer to whether 
MRI is damped or not in the investigated setup, therefore, 
a techinique of separating the stresses should be developed. 
Before attempting that, let us first analyse the behavior of 
the total stresses R^y and M^y measured in our calcula- 
tions. 

The Reynolds stresses of the MRI-unstable calculations 
can well be compared to the forced hydrodynamic ones of 
Sect. 14.2.11 this is done in Fig. |6] where we plot both the 
magnetohydrodynamic stresses (solid lines with symbols) 
and the hydrodynamic stresses (dotted lines) normalised to 
the pure MRI stress. From this figure it can be seen that 
for the sets fcf/fci = 5 and 10, the weakest forcing ampli- 
tudes produce Reynolds stresses that are smaller in mag- 
nitude than the pure MRI-generated stress; this indicates 
that the MRI becomes reduced when external forcing is 
added. It is also evident that for the weakest forcings, ex- 
tra stresses in addition to the hydrodynamic runs are gen- 
erated by the action of MRI (the solid lines run above the 
dotted lines). As the strength of the turbulence grows with 
increasing forcing, the magnetohydrodynamic stress grows 
much more slowly than the hydrodynamic one, so that the 
two curves start approaching each other This could be in- 
terpreted as follows: the action of the forcing on the MRI- 
dynamics is to decrease its contribution to the to- 

tal Reynolds stress, while Rxy continues its growth. At the 
strongest forcings and the highest levels of turbulence, the 
hydrodynamic stress is always larger than the magnetohy- 
drodynamic one, but the shape of the curves is almost iden- 
tical; this could be interpreted as the MRI-contribution be- 
ing nearly zero while Rxy continues its growth. These re- 
sults are consistent with the scenario that the MRI becomes 
damped as forcing is increased, but as the curves signif- 
icantly deviate from each other at high turbulence levels, 
this method cannot be used to definitely determine at which 
point the damping occurs. The stresses from hydrodynamic 
runs are consistently higher than those from the runs with 
MRI and forcing. A possible explanation is that the back- 
reaction due to the Lorentz force quenches the Reynolds 
stress in the latter runs. A weak /c-dependence can be ob- 
served: the larger the forcing wavenumber is, the smaller is 
the Reynolds stress which can be interpreted as an effect of 
the decreasing influence of shear and rotation on the system. 

In contrast to the pure MRI case, in all the forced runs 
large-scale dynamo action develops; such a shear dynamo 
has been previously reported in rotating and shearing flows 
e.g. by Yousef et al. (2008). Due to this a mean mag- 
netic field is generated, showing quite irregular behaviour 
in time disappearing and reappearing without any cyclic- 
ity. The strength of the mean magnetic field is maximally 
around 40 percent of the total magnetic field strength in the 
set kf/ki — 10, and minimally of the order of 10 percent 
with the highest forcings. Such a mean field gives a signif- 
icant additional contribution to the Reynolds and Maxwell 
stresses. We remove this contribution using 

{bxby)^(M^-°"''^)^(Bx)(By), (27) 



where the overbars denote horizontal averaging. 

For the Maxwell stress no such reference point exists 
as for the Reynolds stress. In Fig. |7] we simply plot the 
stress normalised to the pure MRI case. The behaviour of 
the stress is quite similar for the three smallest wavenum- 
ber sets: first the stress is growing as the forcing ampli- 
tude and turbulence level is increasing. At around uum, 
the growth saturates. For the highest wavenumber set, the 
stress is at first strongly decreasing, but saturation is again 
reached around uum- For the largest wavenumber set, the 
total Maxwell stress is always below the pure MRI case; 
also for the set k{/ki = 5 the stress is weaker than the pure 
MRI stress for the weakest forcings. This is a clear indica- 
tion of the MRI-contribution to the Maxwell stress being re- 
duced when forcing is included. If it is further assumed that 
the forcing-induced Maxwell stress grows monotonically as 
function of forcing, then the combined effect with a mono- 
tonically decreasing MRI-contribution at higher forcing am- 
plitudes could result in the saturation of the stress; this, how- 
ever, is only speculative, as no reference point is yet avail- 
able. A stronger fc-dependence compared to the Reynolds 
stress can be seen in the Maxwell stress: the larger the 
wavenumber, the smaller is the stress. There is roughly an 
order of magnitude difference between the Maxwell stress 
of kf/ki = 1.5 and 10. This could be due to the de- 
creasing value of the effective magnetic Reynolds number, 
Rm = Urms/iriki). 

The energy contained in the shear flow is pumped into 
the magnetic energy reservoir by the Maxwell stresses; at 
the same time the magnetic field is being sustained by the 
energy from the turbulent forcing via small-scale dynamo 
action. As long as the MRI is capable of pumping the en- 
ergy from the shear flow, the Maxwell stress normalised 
with the square of the rms magnetic field strength should re- 
main more or less constant. If the MRI-activity was reduced 
due to the external forcing, indications of which were pre- 
sented in the previous paragraph, bi-^s would be expected to 
continue growing while the Maxwell stress gets suppressed; 
in Fig. [8] we plot the Maxwell stress normalised with b^^^ 
as function of rms-velocity. For the sets fcf/fci = [1, 3, 5] 
the normalised Maxwell stress is indeed constant for turbu- 
lence levels weaker than the LSA-predicted limiting value, 
but starts decreasing strongly for stronger turbulence. For 
the set fcf/fci = 10 the Maxwell stress shows a monotoni- 
cally decreasing trend, as even the weakest forcing cases are 
very close to u\i„-^. 

Keeping in mind the behavior of the individual stresses 
(Reynolds stress monotonically growing while Maxwell 
stress saturates as the forcing level is increased), one would 
expect the ratio of these stresses, ~M/R, to start strongly 
decreasing as the level of turbulence exceeds u\i^. As can 
be seen from Fig.|9] where we plot the stress ratio observed 
in each set of our calculations, the decrease of the stress ra- 
tio starts somewhat before Mum, and occurs at the same ex- 
ponential rate for all the wavenumbers investigated. For the 
set fcf/fci = 10, for which the widest range of Mmis/uum is 
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Fig. 8 Maxwell stresses --Mxy generated in the MRI- 
unstable runs with varying forcing wavenumber and forc- 
ing amplitude, normalised to 6^j^s of each run. Colors and 
symbols as in Fig.|6l 
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Fig. 9 Maxwell to Reynolds stress ratio —Mxy/Rxy as 
function of Wrms from MRI-unstable runs with varying forc- 
ing wavenumber and ampUtude. Colors and symbols as in 
Fig.lSl 

available, an order of magnitude decrease is observed when 
the weakest and strongest forcings are compared. Compar- 
ing these findings to the linear results with ordinary Ohmic 
diffusion rj (Sect. 14.1b . indicating that the stress ratio de- 
creases as T] increases, this behavior is indicative of MRI 
suppression. 

To summarise, comparing the stresses developing in the 
forced, rotating and shearing runs to the pure MRI case does 
not directly give information on the suppression or enhance- 
ment of the MRI, as also the external forcing gives rise to 
non-zero stresses. A more appropriate way of interpreting 
the stresses is either to normalise them with the rms-values 
realised in the calculations, or to study the behavior of the 
stress ratio. With the help of all these diagnostics, it seems 
evident that the angular momentum transport due to MRI 
is reduced by the turbulence from external forcing, and that 
the functional dependence on the wavenumber seems to fol- 
low the LSA-FOSA prediction. To definitely conclude that 
the MRI is totally damped, a method of separating the for- 
cng versus MRJ-generated stresses should be found. This is 
tried in the next subsection. 



4.2.3 MRI-stable runs 

We attempt to separate the stresses generated by the ex- 
ternal forcing from the MRJ-generated ones by repeating 
the calculations with MRI-stable shear profile for which the 
angular velocity is increasing outwards (the sign of the q- 
parameter is reversed, so that q ~ —3/2). All the other 
parameters are kept unchanged; the results from this set of 
calculations are summarised in Table |4] 

As in the MRI-unstable calculations, the MRI-stable 
calculations show large-scale dynamo action due to which a 
mean magnetic field is generated. The negative-q dynamo 
appears to be more efficient than the positive-q one; the 
strength of the mean magnetic field reaches 60 percent of 
the total magnetic field strength in some cases. To calculate 
the turbulent stresses, we have removed the contribution of 
the mean field using Eq. dZTl l. 

For the weakest forcings, the Reynolds stress is small 
and negative, showing similarity to the hydrodynamic pos- 
itive q runs presented in Sect. 14.2.11 For stronger forcings, 
positive values, monotonically increasing with the forcing 
amplitude, are obtained. The Maxwell stress is always neg- 
ative, and monotonically increasing in magnitude, but be- 
ginning to show sings of saturation for the largest forcing 
amplitudes. We use these stresses to normalise the ones of 
the MRI-unstable runs; the results are plotted in Figs. [TO] 
and[TT] The trend is quite clear for both of the stresses: for 
the weakest levels of turbulence, there is an excess both of 
Reynolds and Maxwell stresses in the MRI-unstable versus 
the MRI-stable case, indicative of ongoing MRI-activity. 
As the forcing is increased, especially the Reynolds stress 
quickly levels off to a roughly constant value, not showing 
clear correlation with the uiim- The saturation level of the 
Reynolds stress is close to one only in the set with the largest 
wavenumber k[/ki = 10, and increasing with decreasing 
k{. This can be attributed to the asymmetry q ^ —q ob- 
served also in the hydrodynamic calculations (Sect. 14.2. 11 1. 
The Maxwell stress also shows a monotonically decreasing 
trend, the level-off occurring roughly at lium. The saturation 
level, again, is close to unity only for the set fcf/fci = 10, 
and is increasing with decreasing wavenumber We relate 
this behavior to the q ^ —q asymmetry; due to this asym- 
metry, no definite conclusion about the damping limit of the 
MRI can be drawn even with the help of these reference 
calculations. 

5 Application to the Milky Way 

In addition to accretion disks, the disks of spiral galaxies are 
also subject to the MRI, as the angular velocity is decreas- 
ing outwards so that the rotational velocity is roughly con- 
stant as a function of radius outside the very centre (called 
as flat rotation curve, for which q — 1). Sellwood & Bal- 
bus (1999) proposed that the MRI could be responsible of 
the anomalous velocity dispersions observed in the outer re- 
gions of some galaxies, and later, in a series of papers, Pio- 
ntek & Ostriker showed that the MRI can really work in the 
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Table 4 Nonhelically forced 64"^ runs without MRI (shear profiles with angular velocity increasing outwards. Srms stands 



for the total rms magnetic field, 
the mean fields removed. 



{B) the mean magnetic field, and (bxby) are the turbulent stresses with the contribution of 
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Fig. 10 Reynolds stresses from the MRI-unstable runs 
normalised by stresses generated in the MRI-stable runs 
with corresponding forcing wavenumber and amplitude. 
Colors and symbols as in Fig.|6] 



Fig. 11 Maxwell stresses from the MRI-unstable runs 
normalised by stresses generated in the MRI-stable runs 
with corresponding forcing wavenumber and amplitude. 
Colors and symbols as in Fig.|6] 



cool cloudy and clumpy medium subject to thermal instabil- 
ity (see e.g. Piontek & Ostriker 2007). The energy balance 
estimates from the observations of NGC6949 (Beck 2004) 
indicate that magnetic energy could become dominant over 
the kinetic energy in the outer regions of this galaxy, so that 
the energy ratio Em/Ek ~ 3-4. This matches well the en- 
ergy ratios found from numerical simulations of MRI (e.g. 
Liljestrom et al. 2009). 



In the inner star-forming regions of spiral galaxies there 
is another very powerful source of turbulence in addition to 
the MRI, namely supernova explosions (hereafter SN). SN- 
forced turbulent flows have been studied numerically (e.g. 
Korpi et al. 1999; de Avillez & Mac Low (2002); Gressel et 
al. 2008), but no signs of MRI has been reported yet. Here 
we apply the linear analysis result together with the FOSA- 
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prediction to obtain a preliminary estimate whether the MRI 
could be suppressed by SN activity. 

For such an analysis, we need the rotation law for the de- 
termination of il and the expected growth rate of the MRI, 
the SN distribution as a function of radius, and the expected 
rms-velocity and scale of SN driven turbulence as functions 
of the SN forcing. For the rotation law we choose a sim- 
plified fit to the galactic rotation profile similar to the one 
used by Ferriere & Schmitt (2000), shown in Fig. [I2l The 
radial distribution of SNe we obtain from the review article 
of Ferriere (2001), who gives 

ai{R) = (4.8kpc~2Myr"i) 

R — i?0 



xexp 



4.5 kpc 



au{R < 3.7 kpc) 
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x3.55 exp — 



i?- 3.7kpc y 
kpc J 



2.1 



crii(i? > 3.7kpc) = (27kpc"2Myr"i] 



xexp 



(28) 



(29) 



(30) 



(6.8kpc)2^ 

where Rq = 8.5 kpc is the radius of the solar orbit. 

Dib et al. (2006) investigated the velocity dispersion re- 
sulting from SNe of varying frequency in a cube of non- 
stratified and nonmagnetized gas. No rotation or differen- 
tial rotation was included, so the system can be regarded as 
isotropic. SN frequencies of 0.01 to 10 times the frequency 
at the solar radius {Rq — 8.5 kpc) were investigated and 
velocity dispersions, but no information of the scale of the 
turbulence, were reported. From their plots, however, it is 
evident that the scale is approximately 100 pc with the SN 
frequency at solar neighborhood. Gressel (2009) also arrives 
at a similar value by computing velocity structure functions 
from more sophisticated simulations of SN forced stratified 
flows. For our simple analysis, we therefore adopt the value 
l{ = 100 pc, but neglect any variation as a function of SN 
frequency. We then perform an approximate fit to the ve- 
locity dispersion as a function of SN frequency from the 
Figure 13 of Dib et al. (2006), and map the used SN fre- 
quencies over the galactic radius using the SN-distributions 
Eqs. (|28]l - dSOll. The resulting ilm from Eq. (O is plotted 
with a solid line in Fig. [12] 

Finally, from the linear stability analysis, we choose 
A = 1.4, corresponding to the limit where we found all the 
MRI-modes to be damped for q = 1 (see Fig.|2|. We show 
the corresponding damping curve A Fmax = A/2 il{R) 
(see Eq. ( [TS] )) in Fig. [12] with a dashed line. MRI is pre- 
dicted to be damped by SNe in the regions where the solid 
line is above the dashed line. This occurs up to the radius of 
roughly 13.6 kpc. 

The analysis presented above can be improved by in- 
cluding a more realistic description of the turbulent scale 
and rms-velocity and their dependences, e.g. on the inter- 
stellar magnetic field strength and surrounding density, as- 
pects that were totally neglected here. This would require 
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Fig. 12 Galactic application. An approximate fit to the 
galactic rotation profile (dotted line), the resulting Jim from 
Eq. ST4[ (solid line), and the damping curve A Fmax = 
A/2 n{R) (see Eq\l^ (dashed line). 



high-resolution numerical simulations of SN forced flows, 
involving a huge parameter space to explore. Intuitively one 
would expect, however, that the presence of a large-scale 
magnetic field and larger surrounding density would proba- 
bly both result in smaller values of rms-velocities and sizes 
of the individual SN shells. Therefore the ratio Umis/k, ap- 
pearing in the Ohmic diffusion rate in the isotropic case 
(Eq. ([T4]l). could remain unchanged. 



6 Conclusions 

In this study we have investigated the effect of turbulent 
Ohmic diffusion on the excitation and dynamics of the MRI. 
First we presented a linear stability analysis for different 
types of rotation profiles of the form ft cx and pro- 

ceeded by making a prediction in the isotropic homoge- 
neous case of the level of turbulence needed to suppress 
MRI (Eq. 

Next we made an attempt, with the help of three- 
dimensional numerical simulations, to investigate whether 
this simple criterion holds for more complex nonlinear sys- 
tems. We began with purely hydrodynamic calculations 
with both positive (angular velocity decreasing outwards) 
and negative (angular velocity increasing outwards) shear 
parameters q. These calculations gave two important results: 
Firstly, non-zero Reynolds stresses arise solely due to the 
interaction of forcing, shear and rotation, in agreement with 
the earlier results of Snellman et al. (2009), but in appar- 
ent disagreement with the study of Workman & Armitage 
(2008). Secondly, there is a strongly fcf-dependent asymme- 
try in the generated stresses in a transition q ~q, so that 
the smaller the effective forcing wavenumber, the larger the 
asymmetry is in the generated Reynold stress. In the MHD 
regime, therefore, non-zero stresses can occur both from the 
forcing and the MRI, and the separation of the two contri- 
bution is very difficult, as no absolute reference point can 
be established due to the q —q asymmetry. 
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In the magnetohydrodynamical regime sets of MRI- 
unstable (positive q) runs were made, and compared to a 
calculation without external forcing (exhibiting only MRI). 
As anticipated, the interpretation of the results without any 
other reference point than the pure MRI case, was difficult. 
Utilising the corresponding hydrodynamic calculations as 
reference, the MRI versus the forcing contributions to the 
Reynolds stress could be investigated in more reliable man- 
ner From this analysis, indications of MRI damping were 
found. For the Maxwell stress, a normalisation to the rms 
magnetic field strength of each calculation was found to 
be a useful diagnostic, giving further support to the damp- 
ing scenario. The Maxwell to Reynolds stress ratio was 
also observed to be a sensitive diagnostic of the MRI effi- 
ciency. Negative-g calculations were also made in the MHD 
regime; due to the g — g-asymmetry of the generated 
stresses, these runs, however, serve as a rather poor refer- 
ence point. It is possible to study this issue further with the 
help of nonUnear closure models, such as Ogilvie (2003) or 
Snellman et al. (2009). The closure model of the latter study 
could be used to investigate the q — q'-asymmetry in more 
detail, whilst the former could serve as a tool to separate the 
two stress contributions; this, however, is out of the scope 
of this paper. 

All in all, the effect of the external forcing on the exci- 
tation and dynamics of the MRI is to reduce its efficiency as 
angular momentum transport mechanism (seen both in the 
Reynolds and Maxwell stresses, and most clearly in their ra- 
tio), and as the turbulence level exceeds uum, derived by lin- 
ear stability analysis and the first-order smoothing approxi- 
mation, to suppress the MRI-activity. To determine the ex- 
act level of turbulence at which the MRI becomes damped 
in the three-dimensional calculations is not possible due to 
the q — (j-asymmetry, but according to our results, the 
LSA-FOSA-prediction is quite a reliable indicator. 

As evident, the conclusions of our study are in complete 
disagreement with the very similar study of Workman & Ar- 
mitage (2008). The discrepancy is mostly due to the results 
of the hydrodynamical calculations of both studies: our cal- 
culations clearly show non-zero Rl-y, whereas Workman & 
Armitage (2008) find very small B.[.y, which is interpreted 
to be consistent with zero. This likely to be explained by 
the fact that their conclusion is based on only one point al- 
though the stress varies non-monotonically as a function of 
ri and S (Snellman et al. 2009). Although the results of the 
magnetohydrodynamical calculations of both studies are in 
perfect agreement, the conclusions become severely altered 
due to the different reference points. 

One particularly interesting case to apply our result is 
that of the star-forming inner parts of galactic disks, where 
the stellar turbulent energy input produces turbulence via 
stellar winds, radiation pressure and supernova explosions. 
On the other hand, the galactic disks are susceptible to the 
MRI with flat rotation profiles of the type q = 1; is the 
stellar energy input enough to suppress the MRI, and over 
which range of radius does this occur? In this study we made 



a crude estimation utilising the LSA-FOSA-prediction to- 
gether with earlier results from numerical investigations 
of SN forced turbulence (Korpi et al. (1999); Dib et al. 
(2006)) and the radial SN distribution reported in the lit- 
terature (Ferriere 2001). Our analysis gives indication that 
supernova activity can suppress the MRI up to the radius of 
roughly 14 kpc. 
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